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Free- Vibration Analysis of Rotating Beams 
by a Variable-Order Finite-Element Method 

Dewey H. Hodges* and Michael J. Rutkowskif 
U,S, Army Research and Technology Laboratories {A VRADCOM), 

Ames Research Center, Moffett Field, Calif, 

The free vibration of rotating beams is analyzed by means of a finite-element method of variable order. This 
method entails displacement functions that are a complete power series of a variable number of terms. The terms 
are arranged so that the generalized coordinates are composed of displacements and slopes at the element ex- 
tremities and, additionally, displacements at certain points within the element. The displacement is assumed to 
be analytic within an element and thus can be approximated to any degree of accuracy desired by a complete 
power series. Numerical results are presented for uniform beams with zero and nonzero hub radii, tapered 
beams, and a nonuniform beam with discontinuities. Since the present method reduces to a conventional beam 
finite-element method for a cubic displacement function, the results are compared and found to be superior to 
the conventional results in terms of accuracy for a given number of degrees of freedom. Indeed, essentially exact 
eigenvalues and eigenvectors are obtained with this technique, which Is far more rapidly convergent than other 
approaches in the literature. 


Introduction 

R otating beams, which have importance in many 
practical applications such as turbine blades, airplane 
propellers, and helicopter rotor blades, have been studied by 
numerous investigators who used a variety of methods. ’ 
The problem of determining free vibration, response, and 
stability characteristics of nonuniform, pretwisted, rotating 
beams is too complex to solve exactly, especially when 
flapwise and chordwise bending and torsion are considered. 
Thus, most of the various studies found in the literature have 
been concerned with obtaining approximate solutions to 
simplified special cases of the free-vibration problem. It is 
necessary, however, to remember that the methods one uses to 
analyze the simplified free-vibration problems should be 
developed with the more complex problems in mind. Thus, 
such features as reduced computational or programming 
effort, reduced numbers of degrees of freedom, increased 
accuracy, and, above all, ability to incorporate nonlinear, 
nonconservative effects should all be considered important. 

It is also important to identify a set of standard problems 
that can be solved as accurately as possible to provide 
developers of future, more general analyses with a data base 
by which to judge their results. Therefore, the literature was 
examined to see if some standard problems have already been 
identified and solved. 


Discussion of Previous Work 

The investigators examined in the literature considered 
various simplified forms of the rotating beam free-vibration 
problem and obtained solutions for beams with a variety of 
geometries and properties. Among the many methods used to 
solve various simplified forms of the rotating beam free- 
vibration problem arc the Southwell principle,* the in- 
tegrating matrix method, the transfer matrix method,* the 
Mykiestad method,* the Runge-Kutta method,® a 
semianalytic approach based on the Frobenius method,’ two 
modal methods using Legendre polynomials,*’ a mixed 
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variational method, a conventional Rayleigh-Ritz 
method,*'" several forms of the finite-element method, ” 
and an exact solution for some simple configurations.** 
Results presented in these papers were studied carefully with 
the primary purpose of comparison with results from the 
present analysis, which is based on a variable-order finite- 
element method. We will now cite some samples of what was 
found in the literature. 

Lang and Nemat-Nasser used the method of the new 
quotient based on a mixed variational principle to analyze the 
vibration characteristics of a nonuniform blade. This blade, 
which is hinged at the root, was previously studied by 
Wadsworth and Wilde® using Runge-Kutta methods. 
Numerical results presented in Refs. 6 and 10, while in good 
agreement, are not identical, the more significant differences 
being in the moment distribution. Exact numerical results 
obtained from use of the present analysis for eigenvalues and 
eigenvectors, not presented herein, are virtually identical to 
those in Ref. 6. 

Hoa*® has recently presented a finite-element analysis for 
rotating beams with a tip mass. This analysis is based on a 
third-order-polynomial displacement function but is restricted 
to equal-length elements with constant mass and stiffness 
properties. There is, however, as the published comment^ 
indicates, a sign error in the middle term of the elements of 
the centrifugal stiffness matrix in Ref. 16. which leads to 
eigenvalues that are too high. Except for this error, Ref. 16 is 
a very clear and straightforward presentation of a simplified 
finite-clement analysis for rotating beams. The only results 
presented in Refs. 16 and 20 were for uniform beams. 

Unfortunately, in the course of the study it was found that 
several of these references presented a few results which were 
at best of questionable worth. A finite-element analysis based 
on a cubic polynomial-displacement function has also been 
made by Murty and Murthy** for tapered and pretwisted 
rotor blades. When the present finite-clement method is 
reduced to third order, virtually identical eigenvalue results 
(not presented herein) are obtained for various taper ratios at 
zero rotor speed. Unfortunately, this agreement holds only 
for zero rotor speed. All of the results presented in Ref. 13 for 
a rotating beam are in error and for all except the lowest rotor 
speeds the error is substantial. In fact, even the results of a 
Galerkin solution which arc presented by Murty and Murthy 
for comparison are significantly better than the results of their 
finite-element solution. It is obvious that something is amiss 
since the tabulated results attempting to show the convergence 
trend in Ref. 13 do not converge monotonically from above 
for the cases with rotation. 
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This problem is also evident for some of the plotted results 
presented in Ref. 14. The solutions in Refs. 13 and 14 that 
should always be upper bounds are not for some cases. 
Furthermore, whereas the results obtained by using the 
present variable-order finite element will be shown to con- 
verge from above very rapidly to the exact eigensolutions as 
the order of the element and/or the number of elements is 
increased, the same cannot be said for some of the other 
referenced analyses.^ *-’ In Ref. 4, the transfer matrix method 
is used to obtain upper and lower bound solutions for a 
uniform rotating beam. To obtain identical upper and lower 
bounds to seven figures for the first beam natural frequency, 
600 constant-tension element segments are needed. Un- 
fortunately, the converged results with 600 elements are less 
accurate than the integrating matrix results for the same 
beam, also presented in Ref. 4, which more closely agree with 
exact results obtained by using the present analysis. In Refs. 8 
and 9, the displacement is expanded in even* or odd® 
Legendre polynomials. The covergence, for reasons explained 
below, can be relatively slow in such cases. 

In Ref. 17, a simple power series was used to express the 
displacement in each element. Essentially exact frequencies 
and mode shapes were obtained. The resulting mode shapes 
could be differentiated four times and substituted into the 
differential equation and satisfy it to any accuracy desired. 
The numerical results presented in Ref. 17 are for a 
nonuniform beam with discontinuous stiffness and mass 
distributions and are exactly the same as those obtained in the 
present approach for all cases. 

In Ref. 18 the exact solutions for uniform and tapered 
rotating beams were obtained in the form of infinite series, 
and numerical results identical to those obtained by using the 
present method are presented for several configurations. The 
exact solution calculated in Ref. 18 does not apply to general 
nonuniform beams with discontinuities, as in ReL 17 and in 
the present analysis, however. A more complete data base is 
therefore envisioned for the present paper to include essen- 
tially exact solutions for both uniform, tapered, and 
nonuniform beams with discontinuities. 

Present Approach 

It is evident that results presented in the literature fall short 
of a desirable data base in several areas. It is therefore one 
purpose of this paper to tabulate accurate results for a few 
simplified problems, including both frequencies and mode 
shapes for some cases. A second purpose is to present the 
details of derjving the variable-order finite-element for- 
mulation (in particular the shape function) to facilitate its use 
by other investigators. In the present paper, the method of 
Ref. 17 is modified to a true finite-element form so that 
generalized coordinates are actual displacements and slopes at 
various points on the element. This method, of course, yields 
identical results. This variable-order finite-element analysis is 
applied herein to the problem of finding the free-vibration 
frequencies of rotating beams, some of which possess 
nonuniform properties. The method described in this paper 
can also readily be extended to incorporate coupled bending, 
torsion, and extension with geometric nonlinearity and 
nonconservative loading. In Ref. 19, a Galerkin finite-element 
method is applied to a more general nonconservative, 
nonlinear problem. The present method, when applied to that 
problem, will yield equivalent results when the same shape 
functions are used. We note that the shape functions in the 
present analysis, when they are taken to the third order, 
reduce to the conventional cubic shape functions, such as 
found in Refs. 19 and 21 . 

We proceed by first deriving the finite-element equations 
from the principle of virtual work. Then, tabulated results of 
nondimensional eigenvalues as a function of nondimensional 
rotation speed are presented for three representative rotating 
cantilever beam examples: a uniform beam, a uniform beam 
with nonzero hub radius, and a tapered beam. Mode shapes 
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Fig. 1 Geometry of the kth beam element; the angle of the plane In 
which bending occurs, is the same for all elements. 


are also presented for the uniform beam. These examples arc 
chosen as standard problems partly because some of them 
have already been solved exactly. ** In addition, frequencies 
and mode shapes arc also presented for the nonuniform 
discontinuous beam in Ref. 17. Finally, the excellent con- 
vergence obtained by using the present variable-order finite 
element is compared with that of a conventional finite element 
based on a cubic polynomial-displacement function. It should 
be noted that, although all of the results presented in this 
paper are for cantilevered beams, eigenvalues and eigen- 
vectors for beams with other end conditions are also readily 
attainable. ’’ 


Derivation of the Eigenvalue Problem 
Principle of Virtual Work 

We consider a beam element, rotating at constant angular 
speed fl about an axis fixed in space, undergoing bending 
motion described by w in a plane fixed in a reference system 
rotating with the beam (Fig. 1). The beam is assumed to be 
inextensional and the plane in which the beam is bending 
makes an angle S with the angular velocity vector. For ^ = 0 
the motion is purely out of plane (flapping) and for d = r/2 
the motion is purely inplane (lead-lag). The principle of 
virtual work for the Aih beam element is derived in Ref. 17 
and is given by 


f'* /n-, <^^**’* d*v* d6w^ 

Jo V * djr^ * dbr (U 






)cbr=( 


( 1 ) 


where E7* is the distributed bending stiffness, Ti^ the 
distributed tension force, and the distributed mass per unit 
length. The subscript k refers to the Aih element numbered 
from the inboard end to the blade tip. Each clement root is a 
distance at* from the center of rotation and the element length 
coordinate x goes from zero at the root of the element to (j, at 
the tip. The eigenvalue is where = + G^sin^^ and co is 

the natural frequency. The tension force in the Jtth element is 




i: 


(AT* -l-Ar)dAr-f- (d) 


( 2 ) 


where 7\^(f^) = 0 and A/ is the number of elements. Thus, for 
any given beam geometric, mass, and stiffness properties, the 
eigenvalues will be functions of 0. 

Now we define dimensionless parameters 
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( 3 ) 


EL ' 




Elr ' 


^k "^^k^k 

aTz+L Xf^J 


where El, and m, are reference values of bending stiffness 
and mass per unit length and L the total blade length. The 
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principle of virtual work for all elements becomes 

km/ ^ ' 

( 4 ) 


with necessary and sufficient conditions on displacements and 
slopes at element boundaries given as 

w*(^) =w* + ;(0). Sw*(^) =«>*'* + /(0) 



^* + / 


iO), 



h.J 


iO) 


W/(0) —5wj{0)=0 

Wf{0) =6w;(0) =0 (cantilevered case only) 


( 5 ) 


and 



Jn* 


( 6 ) 


where T]^{\)-0. The quantities /??*(»?*), 7*(»?*)» and 
analytic functions and thus the displacement is analytic 
and can be expressed to any accuracy desired by a power series 
in 1 ?*. as illustrated in Ref, 17. Since we are interested in 
developing this method in finite-element form, we will now 
outline the computation of shape functions suitable for this 
purpose. 


where the numbers are the roots of <t>s- 2 iv) = 0 (the Gauss- 
Legendre points associated with the interval). Note that for 
A/= 3 there are no points t;, for 0<i?, <1, Our intent is to alter 
the basis of expansion in Eq, (7) so that the coefficients 
(generalized coordinates) will be w at the points 
y/,yj>y 4 * • • and w' at the points and y^^j. Ob- 
viously, for /^=3 this procedure will yield the standard beam 
finite-element shape functions.^* It is our objective here, 
however, to allow for arbitrary choice of thus enabling 
us to obtain the essentially exact solution as N is increased 
until convergence without refining the clement geometry. 

Now, let 


*v= D (9) 

where ^y(Tt) are the new basis functions. Identification of the 
generalized coordinates at the end of the interval yields 


B/ = w{0) = ^ Aj4>j{0)^ 

J*i 

N4-} N+I 

Bj^w'(0)= £ >4,(-/)^o‘^-y) 

;*/ y=/ 

N-n 

B^ = wU) = Y, A,<t,jU)= T, 

j^i y»/ 

B^^, = W'(1)= ^A^U^-J) (>0) 


Jml 


Shape Functions 

We first expand the displacement w (dropping the bar and k 
temporarily for convenience) in terms of a complete set of 
shifted Legendre polynomials 


s+i 

w= ^ Aj<t>j{rj) O) 

7=/ 

where is the shifted Legendre polynomial of degree 
(/-I) over the interval 0^??^ 1. It should be noted that it is 
imperative for fast convergence that all of the terms be taken, 
not just half as in Refs. 8 and 9. When only half of the terms 
are taken, for example the odd terms, an artificial constraint 
is imposed on the displacement in that the coefficients of 
even-powered terms in an equivalent power series are 
predetermined. It would appear logical to take only the odd 
Legendre polynomials because they are the exact solution to 
the free vibration of a rotating string. The effect of bending 
stiffness on the mode shape of a rotating beam, however, 
cannot be ignored, especially for the cantilever root con- 
dition. A better approach is to take all the polynomials and 
use the boundary conditions to eliminate coefficients of the 
lower degree terms as done in Ref. 17 with a simple power 
series. When this is done, the results obtained are identical to 
those obtained with a simple power series. 

We now identify a set of points in the interval O^y, ^ I 
such that 


y/=y2=o 


For the points in the interval, 

N+/ 

( 11 ) 

A somewhat more convenient expression for programming 
purposes results if we first multiply Eq. (11) by g,4>m(v,) 
sum over /, where is the Gauss-Legendre weight for the 
interval Os TjS 1. 


N-J N-y 

i~i ■=/ 

m=I,2 N-3 (12) 

We now reverse the order of summation on the right-hand 
side and obtain 


N-s ^+/ n-3 

Y^ gi<l>miVi)Bn.2= ^ Aj Y^ (> 3 ) 

,,/ /*/ /-/ 

The inner sum is simply an integral wheny is not loo large 


N-3 


i=I 


s; 




m = 1,2, . . 
m^-j^2N-5 


(14) 


y4^^2 


It is evident that the summation yields the exact value of the 
integral under certain conditions. Thus, 




yN-l^'^N-3 

yN=yN*i = i 


, j^m and js2N— 5 — m 


( 15 ) 
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When y < 2/V- 5 - w, we simply evaluate the sum as it is. The relationship between the coefficients Aj and can now be written in 
matrix form as 

[C]{a]^lD]\b] (16) 

where 


; 

-7 

/ 

. . 



• 

(-/)"' 

0 

2 

-6 

1 

1 

1 

0 1 • 
1 

1 


> 

• 


I 

0 

0 

0 

0 

0 

0 

0 

1/3 

0 

. . . 0 i 

0 

0 

0 

0 

0 

0 

J/5 

. . . o' 

1 


* 






1 

1 

0 

0 

0 

0 





0 

0 

0 

^N-3.NA-f 




1 

0 

0 



0 

0 

0 

1 

0 



/.V+/ 

I 

1 

1 


I 


■ 

7 

0 

2 

6 


p-j 


, 

N{N+I) 




(/V-3) 

10)2X2 ' 


(D] = 

(N-J) x2 

l8j-2^i-2^^j-2 ) 1 (N-3) X (S-3) 

IN-}) x2 

(18) 



(N-i) 

[H2X2 . 



and 


[a] = lA,A,... An^tV = 


S.> = 


S-3 
m = / 


i = N-3,N-2,N-J 
y = N-/,A/,A^+7 


(19) 


(/] = identity matrix [ 0 ] = zero matrix 

Thus, the value of at any point in the interval can be 
obtained through a simple automated process based on 
evaluation of Legendre polynomials at the desired points and 
a single calculation of the above matrices. The integrals 
required in Eq. (4) can be evaluated by Gaussian quadrature. 
If the element property functions 7 * and are polynomials 
in 1 ^, then Gaussian quadrature is exact. It should be noted 
that the Gauss-Legendre points used for the quadrature are 
not the same as those used for the interpolation since more 
points are needed for accurate integration. If C is the 
maximum degree of 7 *, m*, or t*, then the number of Gauss- 
Legendre points Nq required for exact integration is chosen 
such that 


Nc^iG-^D/2 ( 20 ) 


Discretization of the Principle of Virtual Work 
We now substitute the basis functions derived above into 
the principle of virtual work to yield the following eigenvalue 
problem: 


M N+/ N+J -1 

E E*^E«.L(|?w 


with 




B 


NA^l.k _ 

h 


^2,k^I 


^k 

At the root (element I), 




( 2 !) 


( 22 ) 


Bfi —55/; — (? 

B21 = 55^, = 0 (cantilevered root only) 


(23) 


The method of generating basis functions also applies to more 
general problems. If the problem being treated is non- 
conservative, the same basis functions can be used. If it is 
nonlinear, a set of nonlinear algebraic equations would result 
from discretization of the appropriate principle of virtual 
work and more Gauss-Legendre points would be needed for 
evaluating the integrals, but again the same basis functions 
can be used. 


Equation (21) can be expressed as 

hb'^{K-pL^M)b = 0 (24) 

where b is now the column vector of all unknowns 5^* and the 
row vector bb"^ is 55 The eigenvalues are unknown and 
the element stiffness and mass matrices K and M are not yet 
assembled. The global degrees of freedom b are easily ob- 
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Ttbit I Dimensionless frequency ^ ss a function of dimensionless blade angular speed X (cantilever root end, uniform, Xf — 0) 


X 

A/-'5,N=3 

First mode 

Af-1,N=15 

(exact) 

Second mode 

A/=1,N=I5 

M=5,N*=3 (exact) 

A/=5,N=3 

Third mode 

A/=1,N=15 

(exact) 

Q 

3.5160628 

3.5160153 

22.045506 

22.034492 

61.918841 

61.697214 

] 

3.6816930 

3.6816468 

22.191953 

22.181011 

62.062B5S 

61.841763 

2 

4.1373710 

4.1373196 

22.625674 

22.614922 

62.492726 

62.273184 

3 

4.7973618 

4.7972786 

23.330785 

23.320264 

63.202096 

62.984967 

4 

5.5851651 

5.5850015 

24.283711 

24.273349 

64.180850 

63.966760 

5 

6.4498603 

6.4495447 

25.456473 

25.446080 

65.415758 

65.205041 

6 

7.3609370 

7.3603730 

26.819815 

26.809082 

66.891254 

66.683914 

7 

8.3005734 

8.2996369 

28.345571 

28.334083 

68.590234 

68.385951 

8 

9.2582998 

9.2568376 

30.008123 

29.995382 

70.494825 

70.292962 

9 

10.227857 

10.225686 

31.785075 

31.770514 

72.587035 

72.386681 

10 

11.205419 

11.202328 

33.657367 

33.640366 

74.849290 

74.649295 

11 

12.188587 

12.184338 

35.609067 

35.588965 

77.264817 

77.063843 

12 

13.175820 

13.170150 

37.627003 

37.603112 

79.817915 

79.614478 


Table 2 Dimensionless frequency p as a function of dimensionless blade angular speed X (cantilever root end, uniform, X; - I) 


X 

A/=5,N=3 

First mode 

A/=1.N=15 

(exact) 

Second mode 

M=\,N=15 

M = 5.N=3 (exact) 

M=5,N=3 

Third mode 

A/=1,N=15 

(exact) 

0 

3.5160628 

3.5160153 

22.045506 

22.034492 

61.918841 

61.697214 

1 

3.8888709 

3.8888236 

22.385868 

22.375014 

62.263451 

62.043053 

2 

4.8337806 

4.8336888 

23.376585 

23.366042 

63.284581 

63.067548 

3 

6.0820144 

6.0817497 

24.938251 

24.927745 

64.946210 

64.733802 

4 

7.4757038 

7.4750478 

26.968488 

26.957262 

67.194510 

66.986772 

5 

8.9471790 

8.9403581 

29.365944 

29.352835 

69.964967 

69.760710 

6 

10.446342 

10.443866 

32.043670 

32.027244 

73.189437 

72.986335 

7 

11.973160 

11.969071 

34.932917 

34.911582 

76.801539 

76.596449 

8 

13.513664 

13.507389 

37.981697 

37.953793 

80.740275 

80.529532 

9 

15.063167 

15.054077 

41.151554 

41.115408 

84.951845 

84.731533 

10 

16.618936 

16.606363 

44.414262 

44.368224 

89.390169 

89.156329 

1 1 

18.179293 

18.162547 

47.749093 

47.691553 

94.016592 

93.765354 

12 

19.743158 

19.721542 

51.140730 

51.070134 

98.799141 

98.526797 


tained from Eqs. (22) and (23) so that 

^=C6, = (25) 

where C is a matrix that accomplishes the operations in Eqs. 
(22) and (23). Substitution of Eq. (25) into Eq. (24) yields 

6b^C^{K-H^M)Cb = 0 (26) 

Since bb is purely arbitrary, the matrix eigenvalue problem is 
now obtained 

{K~^fi^M)b = 0 (27) 

where the global stiffness and matrices are given by 

K=C^KQ M=C^MC (28) 

For simplicity the assembly process is written as a matrix 
operation in Eq. (28); however, the operations are best 
programmed without matrix multiplication. Since the present 
problem is conservative, both R and A/ are symmetric and M 
is positive definite. The Cholesky decomposition of M into 
LL^ yields an eigenvalue problem involving a single sym- 
metric matrix 

iiUc=L-^KL-'^c (29) 

in which is the eigenvalue and c^L^b. The eigenvalues of 
may be obtained and, if necessary, the eigen- 
vectors may be obtained through appropriate transformation. 
The eigenvalues and eigenvectors may be obtained to any 
accuracy desired by choosing N sufficiently large. Unlike 


the method of Ref. 17, however, the matrices are now well 
conditioned so that double precision arithmetic is un- 
necessary. They are also written in a form convenient for 
finite-element type programming. Unlike conventional finite- 
element methods, however, the displacement at any point on 
the beam can be calculated to sufficient accuracy to achieve a 
smooth representation of the exact waveforms. 

Results and Discussion 

Tables 1-3 present dimensionless modal frequencies ^ at a 
range of dimensionless rotation speeds X of 0-12 for three 
rotating cantilever beams: uniform, uniform with nonzero 
hub radius (off-clamping), and tapered. For the tapered 
beam, a rectangular cross section was chosen where the beam 
depth was assumed to be linear in the dimensionless beam 
axial coordinate r. This choice resulted in a linear variation in 
beam mass and a cubic variation in beam stiffness. In each 
table results for the first three modes are tabulated for two 
different solutions: an approximate solution based on a five- 
element model with a cubic polynomial-displacement function 
(i.e., A/= 5, N= 3) and an essentially exact solution based on a 
single element but with a 15th-order shape function (i.e., 
A/=1,N=15). 

The maximum percentage errors between the approximate 
and the exact solutions in the three modes, respectively, were 
0.43, 0.64, and 0.36 for the uniform case; 0.11, 0.25, and 0.36 
for the uniform with off-clamping case; and 0.013, 0.015, and 
0.29 for the tapered case. It is also interesting to note the 
perfect agreement between results shown in Tables 1 and 2 
and those presented to six significant figures in Ref. 18. 

Numerical results have also been obtained for the finite- 
element solution of the nonuniform, discontinuous rotating 
cantilever beam in Ref. 17. It should be noted that numerical 
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Table 3 Dimensionless frequency /t as a function of dimensionless blade angular speed X (cantilever root end, tapered, Xf = 0) 


X 

Af=5,iV=3 

0 

3.8238530 

1 

3.9866770 

2 

4.4368408 

3 

5.0926865 

4 

5.8787721 

5 

6.7434176 

6 

7.6551936 

7 

8.5957062 

8 

9.5542170 

9 

10.524368 

10 

11.50231! 

11 

12.485652 

12 

13.472865 


First mode 

A/=I,N=15 

(exact) 

3.8237849 

3.9866177 

4.4368017 

5.0926669 

5.8787617 

6.7344986 

7.6551392 

8.5955770 

9.5539580 

10.523905 

11.501549 

12.484474 

13.471130 


Second mode 

Af=I,N=15 

Af=5,N=3 (exact) 


Third mode 

M=^5,N^3 (exact) 


18.325408 

18.317261 

18.481991 

18.474006 

18.944151 

18.936627 

19.690743 

19.683900 

20.691199 

20.685158 

21.910547 

21.905325 

23.313753 

23.309272 

24.868615 

24.864719 

26.547188 

26.543658 

28.326127 

28.322692 

30.186399 

30.182744 

32.112711 

32.108482 

34.092865 

34.087675 


47.403178 

47.264827 

47.554824 

47.417284 

48.006789 

47.871619 

48.750435 

48.619011 

49.772213 

49.645637 

51.054771 

50.933807 

52.578213 

52.463260 

54.321325 

54.212420 

56.262644 

56.159489 

58.381296 

58.283293 

60.657577 

60,563880 

63.073320 

62.982874 

65.612065 

65.523654 


Tapered beam parameters: m = 1 -(/■/2); 7= [1 -(r/2)]^. 


Table 4 Dimensionless frequency /i for dimensionless blade angular speed X = I 
(cantilever root end, nonuniform, Xf = 1/19) 



M 

First mode 

Second mode 

Third mode 

Fourth mode 

3 

2 

1.0691352 

2.6354387 

5.1431835 

14.020328 

4 

2 

1.0663536 

2.5942263 

5.1223762 

9.4830925 

5 

2 

1.0660366 

2.5909878 

4.7319910 

8,8317222 

6 

2 

1.0660098 

2.5909405 

4.7279302 

7.7417328 

7 

2 

1.0660084 

2.5907631 

4.7216555 

7.7388390 

8 

2 

1.0660084 

2.5907040 

4.7189516 

7.6753688 

9 

2 

1.0660084 

2.5906958 

4.7187634 

7.6693676 

10 

2 

1.0660084 

2.5906957 

4.7187630 

7.6692997 

Exact 


1.0660084 

2.5906956 

4.7187589 

7.6691747 

3 

5 

1.0691310 

2.5966530 

4.7313324 

7.7249582 

3 

10 

1.0661534 

2.5909765 

4.7200226 

7.6763372 

Configuration of nonuniform rotating beam with discontinuities: 



r 


rh 



7 

0.05srs0.2 


I 


0.0014«(I+i;)^ 


0.2srs 1 


5[l-(r/2)l 


0.0146(1 +^,)<[1 

-(3r/2) + (3r^/4)l 


values of and EI^ are arbitrary. However, for the purposes 
of the present paper these reference parameters have been 
chosen to give X= 1. In Table 4 the dimensionless frequencies 
for the first four modes arc presented at X = 1 for A/=2 with 
N varying from 3 to 10. As /V is increased the solutions arc 
seen to converge very nicely from above to the exact solution 
that was obtained by taking A/=2 with a sufficient number of 
terms {N-\5) to insure convergence (A/ =5 and N=\0 will 
also produce the exact solution). In addition. Table 4 also 
presents results based on a cubic polynomial-displacement 
function for both a five- and a ten-element model. For A/=5 
and A/=3, the first three modes are in error by less than 
0.30<^^o, and even the fourth mode is in error by only 0.73%. 

In Tables 5 and 6 the displacement, slope, moment, and 
shear are tabulated for the first mode shape of the uniform 
beam (X=l) and the nonuniform discontinuous beam, 
respectively. Here the displacement is normalized for 
deflection of unity at the tip. Results are presented for A/=5 
and and M = 1 and N=\5 in Table 5 and for A/= 5 and 
N= 3 and Af = 2 and N=\5\n Table 6. It is interesting to note 
that the moment and shear are discontinuous for the con- 
ventional results (A/=5 and N=3), as shown in Fig. 2. The 
reason for this is that the cubic displacement field is not 
general enough when differentiated to represent the actual 
variations in these quantities. When yv=3, the natural 
boundary conditions on moment and shear at clement 
boundaries and at the beam tip arc not satisfied. This is the 
reason for the appearance of two numbers under the A/=5 


and N= 3 column for certain values of r. These values of r are 
at the element boundaries and there are thus two values of 
shear and moment (as in Fig. 2) at these points. By simply 
increasing N, however, the appropriate values of moment and 
shear are approached and the results may be made as accurate 
as desired. In fact, the residuals obtained by substitution of 
the first mode shape into the differential equation arc on the 
order of 10 for the uniform beam with X = I for A/ = 1 and 
N=15. Similar results are obtained for higher modes with 
slightly less accuracy. 

The excellent covergence obtained with the present variable- 
order finite element is shown in Fig. 3, where the percentage 
error in the first mode for the uniform rotating beam at X- 10 
is plotted as a function of the number of degrees of freedom 
Mx(N- 1) considered in the analyses. The convergence for 
one element (A/= 1) and two elements (A/= 2) when the order 
N of the elements is increased is plotted along with the con- 
vergence obtained with the third-order element (^=3) when 
the number of elements is increased. The figure clearly shows 
that, for a given number of degrees of freedom, a single high- 
order element is superior to several low-order elements. 
Furthermore, the convergence displayed by the variable-order 
finite element is seen to be considerably faster than that 
obtained with the conventional finite element based on a cubic 
polynomial-displacement function. Similar, although slightly 
slower, convergence is also obtained for the higher modes. 
This convergence is typical of that seen for the other rotating 
beam cases as well as for beams with other end conditions. 
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Table 5 Uniform cantilever beam displacement, slope, moment, and shear distributions for first mode, ir - 1 



r 

w 


dw 

dr 

d^ w 

d d^ H* 
dr^ dr^ 

dw 

dr 

A/-5, 

N-3 

A/-1, 

N=15 

Af-5, 

N=3 

A/-1, 

N=15 

A/-5, 

N=3 

A/=l, 

N=15 


A/=l. 

N=15 

0 0 

0.000000 

0.000000 

0.00000 

0.00000 

3.57027 

3.57602 

-5.14507 

-5.31785 

0 1 

0.016994 

0.017001 

0.33130 

0.33130 

3.05576 

3.05289 

-5.30906 

-5.31007 

0.2 

0.064545 

0.064545 

0.61115 

0.61115 

2.54126 

2.54713 

-5.43842 

-5.25797 






2.53902 


-5.02693 


0.3 

0.137567 

0.137578 

0.84139 

0.84139 

2.06566 

2,06138 

-5.11641 

-5.12358 

0.4 

0.231245 

0.231245 

1.02428 

1.02428 

1.59230 

1.60137 

-5.16378 

-4.87569 






1.58716 


- 4.47171 


0 5 

0.340935 

0.340954 

1.16279 

1.16280 

1.18300 

1.17552 

-4.47756 

- 4.48946 

0.6 

0.462456 

0.462456 

1.26089 

1.26089 

0.77885 

0.79465 

-4.44500 

-3.94607 






0.77156 


-3,28293 


0 7 

0.591922 

0.591953 

1.32364 

1.32366 

0.48362 

0.47158 

-3.21698 

-3.23218 

0.8 

0.726225 

0.726225 

1.35761 

1.35761 

0.19567 

0.22078 

-3.12382 

-2.33920 






0.18701 


- 1.35962 


0.9 

0.862735 

0.862779 

1.37073 

1.37075 

0.07549 

0.05804 

- 1.24547 

-1.26243 

I.O 

1. 00000 

1.00000 

1.37271 

1.37271 

0.03604 

0.00000 

-1.11525 

0.00000 


Table 6 

Nonuniform discontinuous beam displacement, slope, moment, and shear distributions for first mode, y — 1 





dvv 


w 

d / d^ w 

^ r— 


w 

dr 

. O. 1 

Ki 

dr V dr^ 

) dr 


Af=5, 

Af = 2, 

A/=5, 

A/=2. 

Af=5, 


A/=5, 

A/=2, 

r 


N= 15 

N=3 

;v=i5 

N=3 

N-15 

N=3 

15 

0.05 

0.000000 

0.000000 

0.00000 

0.00000 

0.0282493 

0.0499696 

-0.24396 

-1.64836 

O.I 

0.020705 

0.026481 

0.75857 

0.83889 

0.0160511 

0.0096621 

- 1.44882 

- 1.64779 

0.2 

0.123681 

0.124995 

1 .02247 

1.04676 

-0.0083453 

0.0010054 

- 1.85264 

-1.63936 






-0.0039192 


- 1.62697 


0.3 

0.227581 

0.230245 

1.05368 

1.05896 

-0.0023128 

0.0012433 

- 1.55672 

-1.55135 

0.4 

0.334047 

0.336841 

1.07377 

1.07296 

-0.0011034 

0.0010515 

- 1.42773 

-1,41875 






-0.0013018 


- 1.42169 


0.5 

0.442236 

0.444813 

1.08957 

1.08628 

-0.0009224 

0.0008113 

-1.25177 

-1.24695 

0.6 

0.551870 

0.554053 

1.10265 

1 .09828 

-0.0006338 

0.0006082 

-1.04633 

-1.04152 






-0.0006654 


-1. 0461 8 


0.7 

0.662712 

0.664417 

1.11381 

1.10871 

-0.0004635 

0.0004416 

-0.81203 

-0.80815 

0.8 

0.774555 

0.775726 

1.12265 

1.11704 

-0.0003140 

0.0002829 

-0.55513 

-0.55273 






-0.0002974 


-0.55582 


0.9 

0.887116 

0,887714 

1.12790 

1.12205 

-0.0001212 

0.0001096 

-0.28263 

-0.28128 

1.0 

1.00000 

1.00000 

1.12910 

1.12315 

-0.0000300 

0.0000000 

-0.00148 

0.00000 



Fig. 2 Shear force distribution for the first mode shape of the 
uniform cantilever beam, A = 1. 


The trend that clearly emerges is that, for a given number of 
degrees of freedom Mx(A^-l), the method yields more 
accurate results when M is taken as small as possible and N is 
chosen as large as necessary. This trend has been observed in 
other similar approaches to finite-element modeling. 

It should be emphasized that, although the examples 
considered in the present paper have been limited to rotating 
cantilever beams, the analysis method described herein is not 
restricted to these problems. The method is equally applicable 


to rotating and nonrotating beams with other end conditions 
and is readily extendable to include nonlinear coupled bend- 
ing, torsion, and extension with pretwist, transverse shear, 
warp, etc. In fact, this analysis method can also be extended 
to nonconservative problems. Bailey has recently applied 
Hamilton’s law of varying action to Beck’s follower-force 
problem. With the present analysis method, equivalent to that 
of Ref. 23, the exact solution for the critical load of Beck’s 
problem has been obtained to eight places with a single 
element and N- 10. 

Conclusion 

A variable-order finite-element method has been presented 
with application to the free vibration of rotating beams. The 
excellent convergence properties obtainable with this method 
have been demonstrated for several cantilever beam examples 
which include hub offset, taper, and nonuniform, discon- 
tinuous mass and stiffness properties. Where possible, the 
published exact solutions are compared with the present 
results and are found to be identical. In addition, comparison 
has also been made with results obtained by using a con- 
ventional finite element based on a cubic displacement 
function. The results clearly indicate that, for a given level of 
accuracy, as few elements as possible should be used and the 
order of the elements should be large enough to yield the 
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desired accuracy. Moreover, it has been found that with 
respect to computation time, fewer higher-order elements are 
better than more lower-order elements for the same level of 
accuracy. Furthermore, although computer storage 
requirements using the present method have not posed a 
problem for the linear analyses carried out to date this 
question will, of course, have to be addressed for future 
nonlinear analyses. 

The variable-order finite-element method described in this 
paper is currently being extended to incorporate coupled 
bending, torsion, pd extension with geometric nonlinearity 
and nonconservative loading. Preliminary calculations for 
nonlinear and nonconservative problems indicate the same 
favorable convergence properties 
exhibited for the present linear, conservative problem. 
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